---
title: "Untitled"
format: html
editor: visual
---

```{r setup}
#install.packages("nhdplusTools")
library(sf)
library(stars)
library(tidyverse)
library(tmap)
library(rstac)
library(gdalcubes)
library(nhdplusTools)
```

```{r}
#download NHD information for Northern California/Southern Oregon region
download_nhdplushr("/workspaces/cloud-native-geo-avikertesz/NHD", "1801", download_files = TRUE)
```

```{r}
#define file pathway object and view layers in NHD .gdb
GDB_path <- "/workspaces/cloud-native-geo-avikertesz/NHD/18/NHDPLUS_H_1801_HU4_GDB.gdb"

layers <- st_layers(GDB_path)

```

```{r}
#call sub-watershed boundaries from WBD and join to create watershed
wbd <- read_sf(GDB_path, layer = "WBDHU12") |>
  filter(str_detect(HUC12, "^18010106")) |> 
  filter(st_is_valid(Shape))

wbd_outline <- wbd |> st_union()

tmap_mode("view")  
tm_shape(wbd_outline) + tm_polygons(fill = NA)
```

```{r}
#call the flowlines layer
flow_lines <- read_sf(GDB_path, layer = "NHDFlowline")

#join stream order data to flow_lines
huc4_gdb_file_path <- "/workspaces/cloud-native-geo-avikertesz/NHD/18/NHDPLUS_H_1801_HU4_GDB.gdb"

vaa <- read_sf(dsn = GDB_path, layer = "NHDPlusFlowLineVAA")

flowlines <- dplyr::left_join(flow_lines, vaa, by = c("NHDPlusID", "ReachCode", "VPUID"))

#drop M and Z dimensions
flowlines_xy <- st_zm(flowlines)

#clip flowline data to South Fork Eel watershed boundary
Eel_flowlines <- st_intersects(flowlines_xy, wbd_outline, sparse = FALSE)
Eel_flow <- flowlines_xy |> filter(Eel_flowlines)
```

```{r}
tmap_mode("view")  
tm_shape(wbd_outline) + tm_polygons(fill = NA) + tm_shape(Eel_flow) + tm_lines(col = "StreamOrde", lwd = "StreamOrde")
```

```{r}
st_write(Eel_flow, dsn = "eel_flow.fgb", driver = "PMTiles")
```